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FOREWORD 

\ 

»  fhix  report  presents  two  Fortran  rayfracing  axles  anti  related  theory  which  are 
designed  to  treat  underwater  sound  propagation  in  the  presence  ot  range-dependent  sound 
velocity  profiles  tSVPs)  Particular  emphasis  is  placed  on  the  average  horizontal  ray  velocity 
and  its  variation  as  a  ray  traverses  different  SVPs  Efficient  calculation  of  this  variation  is 
of  great  importance  for  improved  modem  multisensor  localization  techniques.  •''The  follow¬ 
ing  pages  discuss  the  adiabatic  invariant  approximation  (AIA)  and  the  use  of  the  codes 
RAN  I  and  RAY2  in  sufficient  depth  that  the  reader  may  apply  them  to  his  own  problems 
or  modify  them  tor  this  purpose  I  his  work  was  supported  by  the  Naval  Electronic 
Systems  Command  PMI  124  under  the  Performance  Evaluation  and  Prediction  program 
at  the  Naval  Ocean  Systems  Center  Portions  of  this  work  have  appeared  in  the  journal 
of  the  Acoustical  Society  of  America 

OBJECTIVE 

I  he  objective  ot  the  work  »as  to  calculate  efficiently  the  variation  of  the  average 
horizontal  ray  velocity  as  the  ray  traverses  different  vnind  velocity  profiles.  Efficient 
valculat  on  of  this  variation  is  ot  great  importance  tor  improved  tnixlem  multisensor 
localization  techniques 

RESULTS 

-  The  adiabatic  invariant  approximation  agrees  well  with  exact  raytracing  results 
even  for  strongly  inhomogeneous  sound  velocity  profiles  Simple  algorithms  arc  developed 
for  demonstrating  this  and  for  computing  the  adiabatic  invariant  in  arbitrary  profiles  for 
a  user-selected  I  amity  of  rays  - 

RECOMMENDATIONS 

A  set  of  curves  or  tables  shoukl  be  prepared  showing  the  functional  relation  between 
the  adiabatic  invariant,  the  average  sound  speed,  and  the  ray  angle,  lor  areas  ot  high 
interest  in  acoustic  surveillance  1  hese  would  provide  a  data  base  for  improved  multisensor 
localization  and  would  reduce  the  calculation  of  average  ray  speeds  to  a  t a ble-look -up  proce¬ 
dure  once  the  location  of  the  source  is  known  approximately 
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I.  INTRODUCTION 


Several  important  problems  in  underwater  sound,  ionospheric  propagation,  and 
seismic  research  involve  the  presence  of  ducts  or  waveguides  which  effectively  trap  the 
energy,  greatly  reducing  transmission  loss  over  the  corresponding  case  of  an  isotropic 
medium  Often,  however,  the  duct  varies  with  range,  precluding  exact  normal  mode  solu¬ 
tions  and  affecting  signal  speeds  in  a  generally  complicated  fashion  The  range  variation 
may  consist  of  changes  in  the  physical  dimensions  of  the  duct  or  in  the  index  of  refraction 
across  the  duct. 

Simple  computational  methods  arc  desired  for  estimating  the  influence  of  such 
range-dependence  The  requirements  basically  are  that  the  procedure  be  readily  checked 
against  known  exact  solutions,  that  it  be  relatively  inexpensive  to  implement  using  archival 
data  bases  and  that  revisions  of  the  data  should  be  readily  incorporated  in  the  method, 
without  requiring  massive  rccompulation 

An  attractive  candidate  for  meeting  these  needs  is  the  adiabatic  invariant  approxi¬ 
mation  <  A! A)  This  presupposes  that  in  the  cases  ol  interest  any  range-dependent  variations 
take  place  over  distances  large  compared  to  a  single  ray  cycle  In  other  words,  the  proper¬ 
ties  of  the  medium  should  be  sensibly  constant  over  a  cycle  Under  such  conditions,  which 
nearly  always  apply  for  deep  sound  channels,  it  turns  out  that  ray  paths  evolve  in  range  so 
as  to  keep  j  certain  quantity  tJl  constant,  known  as  the  action  integral  or  phase  integral  * 
This  "conservation  law"  allows  one  to  follow  a  ray.  in  essence,  without  ever  tracing  it  over 
its  entire  path  It  is  only  nccessars  to  know  in  each  local  region  how  the  quantities  of 
interest  (eg.  maximum  depth  of  the  ray .  vertexing  speed,  average  horizontal  velocity,  etc  I 
arc  related  to  the  adiabatic  invariant  J.  something  easily  done  simply  by  tracing  a  few  ray 
cycles  in  each  local  region 

Although  we  have  been  speaking  in  terms  of  rays,  adiabatic  invariants  appear  in 
mode  theory  as  well  I  ach  mode  corresponds  to  a  particular  vertexing  depth  and  speed  for 
the  associated  family  of  rays  In  referring  only  to  rays,  as  we  shall  do  throughout  whai  fol¬ 
lows.  we  are  assuming  the  civic  t  supports  so  many  modes  that  they  comprise  a  continuum  in 
effect,  their  discrete  nature  may  be  disregarded  without  significant  loss  of  accuracy  Since 
deep  sound  channels  are  of  the  order  of  ten  or  more  wavelengths  across  at  frequencies  of 
interest  for  time  difference  fixing,  the  mode-continuum  approximation  is  generally  quite 
good  Pie  low  order  modes  are  not  strongly  excited  by  acoustic  sources.  which  are  usually 
far  from  the  sound  channel  axis 

When  one  has  only  low  order  modes,  the  adiabatic  invariants  are  the  mode  numbers 
themselves  This  is  intuitively  plausible,  for  if  one  considers  an  arbitrarily  gradual  evolution 
of  one  duct  into  a  second  one.  having,  say.  a  somewhat  different  index  profile,  then  it  is 
clear  that  the  power  initially  in  a  particular  mode  must  equal  the  power  in  that  mode  in  the 
altered  duct  Only  rather  strong  perturbations  will  cause  mode  coupling,  in  their  absence, 
there  is  no  cause  for  energy  to  distribute  itself  among  any  other  discrete  modes 


I  Weston.  IH  .  <»wdrd m  a. Skmh  f  arving  Medium.  Pro,.  Phys  Soc  I l.undnn)v  73. 
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Exceptions  to  this  simple  picture  obviously  exist  in  cases  where  one  duct  evolves 
into  two  In  such  cases  one  cannot  treat  the  splitting  of  acoustic  energy  accurately  unless 
the  details  of  the  SVP  evolution  are  considered  Fortunately ,  such  cases  are  rare  The  proce¬ 
dures  described  in  this  report  should  be  applicable  to  all  ocean  areas  of  surveillance  interest 

The  computer  codes  described  here  sene  two  purposes  RAY  I  allows  one  to  set  up 
a  radically  inhomogeneous  medium  whose  SVP  has  a  prescribed  functional  form  and  to 
determine  numencally  by  how  much  the  approximation  J  *  constant  fails  for  various  hori¬ 
zontal  gradients  RAY  2  allows  one  to  input  an  arbitrary  but  range-independent  SVP  and  to 
trace  j  family  of  rays  in  order  to  find  the  relation  between  J  and  other  properties  of  interest 
lor  rays  in  that  particular  SVP  The  former  code  is  ii  primarily  theoretical  interest,  and 
serves  to  establish  the  accuracy  for  the  AIA  in  analytically  tractable  cases,  while  it  is  our 
hope  that  the  latter  will  serve  as  a  model  tor  implementation  ol  AIA  methods  in  acoustic 
surveillance 

In  section  II  we  describe  in  general  terms  the  scope  of  RAY  I  and  RAY2.  we  review 
briefly  the  adiabatic  invariant  theory  as  applied  to  acoustic  ravs.  the  analogous  problem 
in  classical  mechanics,  and  finally  give  the  AIA  signal  speed  algorithm  Sections  III  and  IV 
provide  instructions  in  detail  on  the  procedure  for  executing  RAY  I  and  RAN  2  In  section  V. 
we  shall  treat  the  internal  logical  structure  ol  RAN  I  and  RAN  2  in  order  to  pave  the  way 
for  any  changes  deemed  necessary  or  desirable  I  istings  ol  the  Fortran  source  code  appear 
in  sect  Kin  VI 


II  GENER  AL  DESCRIPTION  OP  RAY  I  ANDRAY2 

In  this  section  we  shall  describe  the  general  features  of  RAN  I  and  RAY2  and  the 
theory  behind  these  codes  In  sections  III  and  IV  we  shall  present  specific  information 
needed  to  make  runs  with  these  codes  Mere  we  are  concerned  with  the  practical  need  for 
such  codes,  and  the  underlying  concept  of  using  an  adiabatic  invariant  as  a  means  for  follow¬ 
ing  the  evolution  of  rays  in  range  dependent  SVPs 

A  THE  ADIABATIC  INVARIANT  FOR  SOUND  RAYS 

fhe  basic  problem  of  time-difference  fixing  provides  a  context  lor  introducing  the 
adiabatic  invariant  concept  For  simplicity,  suppose  that  acoustic  signals  propagate  with 
the  same  speed  in  the  txean  regardless  of  direction  In  this  isotropic  ocean,  a  pair  of  sensors 
which  receive  a  given  sipr  nay  determine,  based  on  the  difference  of  ar  -val  finies.  a  locus 
of  points  on  which  the  •  c  must  Ik-  1  he  locus  is  a  hyperbola  with  the  x-nsors  as  foci 
If  Rj  and  Rs  represent  dis. arises  to  sensors  I  and  2  respectively.  and  c.Nt  represents  the 
equivalent  delay  length,  then  IRj  -  Rsl  *cAl  determines  the  hyperbola  A  third  sensor, 
or  come  additional  data  from  one  of  the  original  two.  then  suffices  to  locate  the  source 
uniquely  (Ambiguities  arising  in  particular  geometries  arc  not  of  interest  here,  although 
they  may  be  important  in  practice  f 

The  procedure  just  described  works  only  approximately  in  actual  ocean  environments 
For  long-range  detections  on  the  order  of  several  hundred  kilometers  the  signal  speed 
generally  will  vary  with  range  Modern  acoustic  surveillance  systems  arc  capable  of  suffi¬ 
ciently  precise  measurement  that  this  inhomogeneity  poses  a  maior  limitation  on  localization 
Mow  can  one  determine  an  appropriate  average  sound  speed  so  as  to  construct  the  needed 
locus’’ 
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A  natural  solution  is  to  use  a  conserved  quantity  for  sound  rays  in  range-dependent 
environments.  Weston'  has  pointed  out  that  the  quantity 

J  *  f  c"'  sind  da  " ) 

is  conserved  provided  that  changes  in  the  SVP  are  sulficienlly  gradual,  or  adiabatic.  Here 
c  is  the  I  oval  depth -dependent  sound  speed,  and  is  regarded  as  being  range-independent  over 
the  penod  of  integration  (one  ray  cycle  i.  d  is  the  ray  angle  measured  with  respect  to  hori¬ 
zontal.  and  z  is  the  depth  coordinate  (usually  taken  as  positive  downwards)  More  detailed 
discussion  of  the  mvanance  of  J  in  range-dependent  SVPs,  as  well  as  a  short  bibliography  on 
the  topic,  will  be  found  in  reference  2  A  proof  ot  the  adiabatic  invariance  of  the  action 
integral  will  be  found  in  reference  .1  At  this  point  it  may  be  ol  interest  to  consider  the 
mechanical  analogue  of  trapped  sound  rays  and  eq  1 1 1 


I .  Uamical  Mechanical  Analogy 


Ray  theory  is  of  course  an  approximation  to  an  “exact"  wave  theory  of  propagation 
In  the  same  sense,  classical  mechanics  is  an  approximation  to  the  more  precise  quantum 
mechanical  theory  In  both  ..axes  the  smallness  of  the  wavelength  w  ith  respect  to  the 
physical  dimensions  ol  the  medium  or  measuring  apparatus  involved  permits  one  to  pass 
over  to  the  geometrical  optics  limit,  ignoring  small  wavelength-dependent  corrections.  In 
classical  mechanics,  the  motion  ot  a  system  is  obtainable  from  Hamilton’s  principle,  which 
is  a  variational  principle,  stating  that  the  actual  trajectory  followed  is  such  as  to  minimize 


the  integral  /  U4-  d  .  t)  dt  t  Actually  the  statement  of  Hamilton's  principle  permits 


maximization  as  well,  but  we  shjll  not  be  concerned  with  this  in  what  follows.)  Here  L  is 
the  lagrangian.  a  function  of  the  coordinates  and  velocities,  symbolically  denoted  q  and 
respectively  ,  and  possibly  ot  time  l  The  coordinates  and  velocities  arc  functions  of  time 
The  system  moves  from  some  fixed  starling  point  (qif  j  I,  qtt  | ))  to  a  final  point  (qtfs). 
q(|s)j.  also  fixed  Hamilton’s  principle  is  written 


ft  /!-  Ldt  -  0  . 

’l 

where  ft  indicates  an  infinitesimal  sanation  ol  the  path  fquation  ( 2)  states  that  the  system 
follows  a  minimizing  trajectory,  since  any  infinitesimal  variation  about  the  trajectory  must 
produce  no  change  in  the  value  of  the  integral 

In  classical  mechanics,  the  equations  ot  motion  follow  if  one  uses  for  the  1  Jgrangian 
|  »  T  -  V.  where  T  is  the  kinetic  energy  and  V  the  potential  energy  of  the  system  Then  the 
calculus  of  variations  yields,  as  a  theorem  tmm  cq  I  2  I.  Newton  s  equation  ot  motion. 


2  Shockley.  RC.  Ptrtxul  anJ  \nnparartal  Rt*  Spmlt  in  St* -iffy  Hmifr  Drprruirril  SOt AR  Chmnrit, 
I  Acmist  Soc  Am  v  M.p  1 1 71-1177  1 
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V .manorial  principle*  (with  different  choice*  for  L)  also  lead  to  the  equation*  of  electro¬ 
dynamics.  hydrodynamic*,  quantum  mechanic*,  etc.,  in  similar  lash  ion 

ll  we  consider  Fermat's  principle,  the  •similarity  to  Hamilton  *  principle  is  clear 
Fermat's  principle  slate*  that  a  ray  will  follow  that  trajectory  minimizing  the  time  of 
flight  between  a  given  initial  and  filial  point 

6  / j  dt  •  6  fjds/c  -  0 

where  d*  i*  an  infinitesimal  element  of  length  along  the  ray  path  and  c  the  speed  of  propa¬ 
gation  hi  the  medium 

For  concreteness.  »r  restrict  our  discussion  to  a  slab  symmetric  medium  and  let  x 
be  the  horizontal  and  /  the  vertical  coordinate  I  hen  eq  (.h  becomes 

6  /‘c'1  (I  ♦z:>‘:d\  »  0  ,4> 

where  i  a  dx/dx  i*  the  slo|K-  of  the  lay  I  hi*  result  now  is  identical  to  eq  <2l.  Hamilton  » 
principle,  under  the  transformation 

t  -*  x  . 

q-z  . 

Lfq.  q.tl-*(l  ♦  xl  ( 

Thus  range-  plays  the  role  ot  time  in  the  mechanical  problem,  and  the  depth  z  is  analogous 
to  particle  displacement  Hence  z  land  is  analogous  to  velocity  Note  that  l  does  not 
resemble  a  quantity  like  I  \ 

It  is  natural  to  ask  what  equation  involving  the  "acceleration"  z  is  analogous  to 
Newton's  equation  ol  motion  The  answer  provides  insight  into  the  behavior  ol  rays  in 
deep  sound  channels  and  bears  directly  on  the  rav tracing  routine  employed  in  RAVI  and 
RAY’  The  filler  l agrangc  equation,  which  follows  as  a  theorem  from  the  variational 
principle  in  cq  l  2l.  as  mentioned  above,  states  that 

td/dtlftL/dq  -dL/dq  *  0  «» 

If  the  reader  applies  the  transformation  given  m  eq  l  $i  then,  alter  winu-  algebraic  mampula 
turn,  csy  fbl  wilt  be  seen  to  become 

V  *  q|  ♦  }‘jc**  dc/dz  ♦  ill 

a  result  obtained  previously  by  Milder  2  4  I  quation  «7>  is  similar  to  Newton's  equation  ot 
motion,  giving  the  acceleration  av  a  function  ot  displacement,  velocity,  and  time 

RAVI  and  R A Y2  use  a  fourth  order  Runge-kntt  i  '•  heme*  to  integrate  cq  C).  and 
thus  obtain  the  ray  depth  as  a  function  ot  range  fquation  <  i  describes  the  ray 

4  Milder  DM.  Rav  ar*/  KVn  r  Imanantt  for  SOUK  Channel  1  Acoust  Sck  Am  .v  4«v 
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curvature  in  terms  of  the  local  sound  speed  denvatives  in  depth  and  range  and  the  ray 
slope  l  or  prescribed  initial  conditions  /<0)  and  /(Ol.  and  a  given  function  c(/.\),  it  speci¬ 
fies  a  unique  trajectory  satisfying  Fermat's  principle 

When  the  medium  is  range-independent,  the  second  term  on  the  right-hand  side  of 
cq  ( 7>  vanishes  A  turthei  simplification  comes  about  since  in  range-independent  media 
Snell  s  law  holds  m  the  form  cv  =  c/c os®,  where  c%  is  the  constant  vertexmg  speed  (where 
/  3  0 1  Thus  eq  ( 7 1  becomes 

7  ■  -  c“(dc/dr)  c-*  3  q3/Ji)l-cj/2c*)  <R) 

Hence  m  range  independent  SVPs  the  motion  of  a  ray  is  precisely  analogous  to  the  trajectory 
of  a  mechanical  particle  of  unit  mass  under  the  action  of  a  potential  field  given  by 


I  his  result  allows  considerable  insight  into  the  trajectories  ol  rays  l  or  example,  if 
we  consider  the  HirsJi  profile 

V  •*  V, 

cU)  ■  cq/(  I  -  t~  *~ I  . 
then  the  potential  is 

A  1  > 

Vt/t  *  -  tcv  cq»*  ( I  -  /-  'a*  l 
■  (c  /cq  al*  /*  ♦  const  . 

which  is  a  quadratic  function  ol  /.  the  same  as  a  simple  harmonic  oscillator.  Ray  paths 
must  therefore  be  sme  waves  in  this  SVP  It  is  simple  to  sketch  potential  functions  tor 
an>  reasonable  SVP.  and  range-dependence  may  be  included  m  this  analogy  by  allowing 
V(/>  to  vary  slowly  with  time 

When  the  medium  is  range-dependent,  then  the  analogous  classical  mechanics 
problem  involves  a  time-dependent  potential,  as  one  tan  sec  from  the  transformation  eq  (5  I 
Hence  the  energy  of  the  analogous  particle  is  not  conserved,  nor  arc  the  positions  of  the 
turning  points  where  q  3  0  or  /  3  0.  If  the  variations  are  gradual,  however,  the  action 
integral 


is  approximately  convened  *  Here  p  =  dL'rfq  .  p  is  tailed  the  momentum  canonically 
conjugate  to  the  coordinate  q  in  classical  mechanics  If  the  reader  again  applies  eq  (5l  and 
uses  /  =  tan#,  he  will  find 

p  3  z/c(l  3  sinfl/c  . 

and  thus  cq  ( I  f  follows  as  the  adiabatic  invariant  for  sound  rays  Note  that  the  integral 
coven  one  cycle  of  the  trajectory  and  represents  the  area  of  an  (a  most!  closed  loop  traced 
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out  in  pn  space,  known  as  phase  space  Ihc  path  is  a  closed  h>op  for  range-independent 
SVP*. 

2-  AIA  Algorithm  for  Average  Horizontal  Ray  Velocity 

One  ma>  use  the  invariance  of )  to  compute  the  average  hon/ontal  velocity  of  a  ray 
traversing  a  range-dependent  SVP  as  follows  We  assume  it  is  possible  to  break  up  the  entire 
region  of  interest  into  local  regions  bounded  by  some  set  of  contiguous  contours  inside  each 
of  which  ctz)  is  a  known  function,  which  cjn  be  regarded  js  range-independent  within  the 
kval  region 

l  irst  one  constructs  tables  showing  the  relationship  between  J  and  c.  the  average 
horizontal  ray  velocity,  for  each  local  region.  Mils  may  require  tracing  perhaps  10  to  50 
rays  tor  each  SVP  fach  ray  need  only  be  traced  lor  a  small  number  of  cycles  (one  or  twol 
to  find  J  and  c  for  the  given  initial  conditions  these  tables  then  allow  one  to  find  how  c 
varies  trom  region  to  region  lor  a  particular  J  value,  by  employing  simple  table-look -up 
procedures 

t  his  method  does  not  yield  detailed  information  on  the  trajectory,  such  as  positions 
of  ray  turning  points  in  range,  but  this  is  the  cost  of  obtaining  estimates  of  c  lor  long-range 
paths  with  simple  algorithms  I  or  a  prescribed  path,  one  then  nuy  find  the  effective  signal 
speed  lor  given  initial  conditions  (depth  and  ray  anglel  from  the  c  values  in  each  region 
according  to  the  fraction  of  time  spent  in  each  region  It  K  denotes  the  total  range.  R,  the 
distance  traveled  in  region  i.  and  c  ,  the  average  horizontal  velocity  in  region  i.  then  (In¬ 
effective  signal  speed  is  given  by 

c  ■  R'S (R.  c  I  (lit 

i  ’  1 

Additional  algorithms  are  necessary  to  select  a  range  of  initial  conditions  that 
characterize  rays  which  carry  most  of  the  acoustic  energy  Depending  on  the  geometry, 
these  may  range  from  simple  inspection  procedures  to  complicated  ray-trace  routines 
requiring  detailed  bathymetry  charts  Reference  2  discusses  possible  methods  lor  selecting 
appropriate  ranges  of  J  values  using  the  output  ol  R.AY'2  (see  section  A  ll  further  remarks 
on  this  are  included  in  section  IK‘.  where  we  indicate  how  certain  of  the  quantities  com¬ 
puted  by  RAY2  facilitate  ray  selection 

B  GENERAL  DESCRIPTION  Of  RAY  I 

A  listing  of  the  code  RAY  I  will  be  found  in  section  \f.  and  instructions  on  its  use-  in 
section  III  Here  wr  shall  describe  the  general  features  ol  RAY  I 

RAY  I  provides  quantitative  answers  to  the  question  of  the  validity  of  the  AIA 
farlier  we  noted  that  J  is  conserved  for  sufficiently  slow  variations  in  the  SVI*  More  pre¬ 
cisely.  J  is  conserved  lor  nearly  sol  provided  that  any  changes  in  the  SVP  arc  fractionally 
small  over  one  cycle  distance  (in  range  I  of  the  ray  This  follows  from  the  fundamental 
law  of  ray  theory.  Fermat's  principle  of  least  time  ’  as  discussed  in  section  II  A  I  We  refer 
the  reader  desiring  further  discussion  of  this  to  landau  and  l  ifshitz-*  and  the  bibliography 
of  reference  2 
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In  real  ocean  environments,  the  jdiabatic  condition  is  satisfied  tor  KR  (Refracted 
Retracted)  rays,  but  one  still  may  ask  what  degree  of  inhomogeneity  is  required  betore  the 
AlA  fails  badly  RAY1  provides  specific  answers  to  this  by  allowing  the  user  to  control  the 
horizontal  gradients  ot  parameters  in  the  SVP  function  Thus  both  the  degree  ol  inhomo- 
geneity  and  the  effect  ot  the  functional  form  c(z  )  may  be  investigated 

Presently  three  self-similar  SVPs  are  available  js  subroutines  in  RAN  1  I  "hey  are 

the  parabolic  profile. 


ctz.  x)  “  Cq(/»  |  I  ♦  /*  a*  < x )| 
the  Ilirsch  profile, 

ctz.  \  )  *  Cgtzl  |  I  -  Z“  a*  tx  1| 
and  the  bilinear  profile. 


(12a) 


(12b) 


ctz.  v  )  *  Cq(  v.  1 1  I  •  z  atx  )1  z  >  0  ( <  0 ) 


1 1  2c  i 


I  igure  I  shows  these  SVPs  tor  range-independent  S\  Ps  I  ach  <>l  equations  (3)  describes  a 
slab  symmetric  (depending  only  on  range  v  and  depth  zi  sound  channel  with  an  axis  at 
z  0.  and  an  axial  sound  speed  c()lxl  Both  c^ivi  Jtul  at  x  i  are  linear  (unctions  of  x  in  the 
current  version  ol  R  \Y1  By  varying  their  derivatives  with  appropriate  input  data,  the  user 
may  study  the  degree  to  which  nonadubjtic  changes  m  the  SVP  affect  the  action  integral, 
or  the  discrepancy  between  the  computed  average  horizontal  ray  velocity  and  the  predicted 
value,  assuming  dj  d\  0  Details  on  the  program  inputs  will  be  deferred  until  section  III 

R  AN  I  allows  the  user  to  select  one  ot  the  three  SVPs  and  control  tracing  of  rays 
in  the  presence  ol  constricting  ( it  da  dx  ■  0)  or  widening  (da  Ox  >  (M  sound  channels 
('•radients  m  c()  and  a  are  independent  In  addition,  the  ray  tracing  routine  employed  (a 
fourth  order  Runge-kuttj  algorithm  applied  to  the  I  uler  I  agrange  equation  ol  motion 
obtained  trv>m  Fermat's  principle)  has  an  automatic  test  and  redefinition  routine  which 
refines  the  step  size  dx  when  there  are  fewer  than  150  steps  per  half-cycle,  although  this 
number  can  be  easily  modified  (see  section  \  i 

for  completeness,  we  present  in  table  I  the  analytic  formulae  tor  rjy  paths  and  the 
average  horizontal  rjy  velocity  in  each  ol  the  SVPs  in  RAN  I  tor  range-independent  con¬ 
ditions  for  the  parabolic  SVP.  ui  and  ftui  i  arc  the  normal  elliptic  integrals  of  the  first 
and  second  kind,  respectively,  with  argument  z  a  modulus  squared  lo-l  I-  2<o--»l )  where 
cv  c().  sn  u | .  cn  ii|  and  dn  uj  are  Jacobian  elliptic  functions,  ft k  I.  K(k (  and  llffl-.  k I 
are  complete  elliptic  integrals  ol  the  first,  second,  and  thud  kind  with  k  -  =  (a-  I  (-  t  2a- *2 1 
and  the  parameter  d*  =  <  I-cosOq)  2  *  t)()  is  the  angle  with  respect  to  horizontal  at  which 
the  ray  crosses  the  sound  channel  axis  at  z  *  0  in  all  formulae 

I  igures  2  and  }  show  the  agreement  found  for  fl()(0i  *  AO  in  an  earlier  study' 
between  predicted  sound  speeds,  assuming  c)J  r)\  0.  and  those  computed  by  RAN  I  in  the 


*  Bvrif.  Pf  ami  Friedman.  Ml).  Hjr*Jh-s>k  of  Uliptif  Intcgralt  ft*  f  nrmrm  and  Sctminlt,  Springer 
Verlag.  New  York  ( l*J?| ) 


presence  of  unoui  gradients  in  atx),  the  depth  scaling  function  Hie  agreement  generally 
is  good  even  when  the  environment  changes  rapidly.  tor  example  by  40  tin  atUI  in  a  single 
ray  vyde  Such  agreement,  althougli  it  vannot  be  expected  tor  all  SVP  choices,  lends  con¬ 
siderable  support  to  the  general  vahdil>  of  the  AlA  t.hiantitative  answers  to  this  question 
are  of  interest  since  the  AlA  method  seems  to  offer  a  simple  solution  to  some  time 
difference  fixing  (  1 1)1  >  problems  data  bases  lor  implementing  AlA  sound  speed  algorithms 
are  easily  generated  the  table-look-up  \|-\  computations  typically  are  lar  less  time 
consuming  than  long  range  ra>  tracing  and  the  method  permits  very  tlexible  updating  ot  data 
in  individual  SVP  regions  without  requiring  long-range  ray  tracing  computations  to  be  repeated 

In  figures  2  and  3  Xp  is  the  cycle  distance  and  x^  the  c ritual  range  at  which  the 
channel  constricts  to  an  extent  such  that  rays  are  reflected  and  vannot  propagate  further 
Hie  notation  xk  x.,  t  x  0 1  is  the  ratio  ot  t fie  .  ritu  a)  range  to  a  v  >  v  le  distance,  ev  alualed  at 
the  initial  range  thus  tor  example,  it  xv  Xp  V  the  critical  constriction  occurs  only  3 
cycle  dictates  from  the  ray  s  starting  point  Ihe  vertical  axis  is  the  normalized  average 
horizontal  ray  speed  tc  *  I  lot  a  ray  on  the  channel  axis!  and  the  horizontal  axis  is  the 
normalized  vhannel  width  parameter  la*  a  a(j  where  a  is  the  s^ale  length  in  the  table  I 
formulas,  and  iq  >'  fhe  initial  value  ot  ai 
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u>  PARABOLIC  SVP 


It '  Bil  ISC  AH  SVP 


fifurr  I 


Thrc*  jrulvticjl  formi  for  modfl  «•  *»m«>  i*)<xitv  profilr*  Jfiom  tq  f  I  m  dimen«onleM  form 


SVP 

function 

c<i> 


Kav  I  tajo.  imv  anJ  Avenge  Speed 


Cyll  ♦/-»*» 


v  a  (2a*  ♦  2)*  "  |(2a*  »  2)  1 4 u 1 1  -  (a*  ♦  2o»Uj 

-  (a  - 1 1*  vi  u  |  in  u  |  dn  u  1 1 

C  vy  It^l.'llU  -  k|k»  |(l  *  ^v«0tll(^:  n 

-  k(ki| 


^  1  u 
H/i-i-vr  • 


/  “  J  UI10Q  O’M  V  J  COtOy) 
C  Cy  2  CIM0Q  I  I  ♦  COI*#q) 


■•yl  I  •  /  1» 
/  <  0(  •  0» 


1/  j  ♦  I  r  ♦  I  v  «  -  unflyl* 


Cy  •  tanfly  lri  *  tanOy) 


Input  ti>  RAN  I  consists  ot  j  control  or  vwiKh  varuMc  indicating  which  SVP  is  to  he 
used.  a  \ct  ot  parameters  del  mini;  the  medium  le  g  axial  sou  ml  speed,  vortical  scale  factor, 
anil  gradients  with  res|H-ct  to  range i  a  set  ot  initial  tav  conditions  (depth  and  slope!,  and 
parameters  lor  the  range  step  si/c  and  maximum  range  Output  is  m  the  form  of  j  summit) 
table  listing  the  fas  coordinates  at  each  sertexing  point  Ihese  are  also  called  lay  turning 
points  f  or  RR  ra\s.  all  turning  points  are  charaetcn/cd  h>  d/  dx  *  (I  (In  RAN  2  the  surtace 
reflection  lor  RSR  (Refracled-Surlacc-Rcflcitcdl  ra>'  leads  to  a  discontinuous  slope  i 
In  addition  to  the  coordinates.  R  AN  I  also  lists  the  maximum  depth  reached,  the  cycle 
distance,  the  phase  integral  evaluated  oser  each  halKyvIc  and  the  average  horizontal  lav 
velocity  I  samples  are  presented  in  scetion  VI 

Although  R  AN  I  has  established  that  the  AIA  works  well  even  in  severely  range- 
dependent  SVPs  ot  certain  kinds  one  should  note  that  only  linear  variations  ot  the  depth 
sealing  factor  have  been  programmed  An  extremely  interesting  question,  which  eould  be 
answered  by  minor  modifications  in  R  AN  I  involves  the  conservation  ot  J  in  two  range- 
independent  SVPs  whieh  are  connected  by  some  inhomogeneous  region  wherr  atx)  vanes 
smoothly  In  each  range-independent  region  the  phase  spase  traiec  tones  are  closed,  and 
J  is  well  defined  as  the  area  enclosed  bv  the  phase  space  orbit  As  R  AN  I  is  currently  written, 
the  trajectories  are  never  closed  (unless  the  user  inputs  da  dx  ■  0  as  a  condition)  A  second 
question  is  the  effect  of  nonsymmetnc  SVPs  AA  ith  minor  modifications.  R  AN  1  could  be 
used  to  study  other  more  realistic  SVPs  than  those  presently  programmed 

C  CifcNf  RAL  DESCRIPTION  OP  RAY2 

RAN  2  consists  of  four  fortran  routines  which  trace  a  family  ot  rays  in  a  user  supplied 
SVP.  computing  lor  each  ray  the  turning  point  coordinates  (depth  range,  and  time),  the 


AIA  PREDICTION 


EXACT  NUMERICAL  RESULTS 
SYMBOL  I  Xt/XplX  •  01 


la  I  PARABOLIC  SVP 


Companion  of  AIA  pi  cJk  tii«m  and  cxa«.l  numetual  itiulli  foj  paiaholk  SVP 
Ifmm  tel  'I 


lb>  MIRSCH  svp 


AIA  PREDICTION 


Fijurf  '  ( nmpanvn  of  AIA  prrdictiom  and  «tki  mimtntal  rrailtt  foi  a  llirach  SVP 
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average  horizontal  ray  velocity  between  turning  points,  the  vertexing  sound  speed,  and  the 
action  integral  computed  over  one-hall  cycle 

In  discussing  the  adiabatic  invariant  J  -  f  (sine? >  d/  c  ( which  is  also  called  the  action 
integral  in  analogy  to  classical  mechanics,  or  the  phase  integral  since  it  represents  an  area  m 
phase  space)  we  outlined  a  general  sc' he  me  lor  efficiently  computing  an  effective  signal  speed 
over  long-range  paths  Essentially,  the  procedure  was  reduced  to  a  table-look-up  method  to 
obtain  from  a  data  base  of  previously  computed  c  vs  J  tables  a  particular  c  value  hi  each 
SVP  region  A  weighted  average  of  these  is  then  computed  to  obtain  an  average  sound 
speed  for  the  particular  J  value  (that  is.  lor  the  selected  initial  conditions  for  the  ray 
trajectory)  over  the  chosen  path 

As  a  lust  step  m  implementing  such  an  algorithm,  see  tested  the  hypothesis  "J  is 
conserved"  for  three  range-dependent  SVPs  of  the  parabolic.  Hirsch.  and  bilinear  types  using 
RAY  I  The  results  tend  to  support  the  assumption  that  J  is  an  invariant  in  range-dependent 
SVPs  to  an  approximation  suitable  for  pmctie.il  signal  speed  calculations  m  typical  ocean 
environments  supporting  HR  rays  * 

RAY’  provides  a  means  lor  calculating  the  relationship  between  c  and  J  lor  arbitrary 
SVPs  of  the  user's  choice  Both  RR  and  RSR  lays  are  allowed  and  typical  output  can  he 
presented  as  in  figure  4  In  addition  to  the  c  vs  J  curves,  one  also  needs  a  means  lor  selecting 
an  interval  of  J  values  characterizing  the  rays  which  carry  most  of  the  acoustic  signal  energy 
Ibis  may  be  done  by  eliminating  rays  which  exceed  a  cut-off  depth  or  by  including  only 
those  rays  launched  from  the  source  (or  arriving  at  the  receiver)  within  a  fixed  angular 
window  RAV2  provides  a  list  ol  lay  turning  point  depths  and  a  summary  of  the  initial 
conditions  to  facilitate  this  selection 
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figure  4  Relation  between  C  and  I  fot  a  selection  of  SVPt  tret 
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Further  discussion  ot  RAY2.  including  input  data  requirements  and  internal  logic, 
will  be  found  in  sections  IV  and  V.  It  is  worthwhile  to  stress  here  that  RAY 2  treats  rays 
only  in  homogeneous  (range  independent!  SVPs.  and  traces  a  family  of  40  rays  covering 
0.5°  to  20°  in  depression  angle  at  the  source  out  to  the  fifth  vertcxing  point  (The  latter 
parameters  are  simply  adjusted  by  internal  code  changes.)  Rays  are  traced  over  several 
hall-cycles,  even  though  a  single  half-cycle  suffices  to  find  c  and  J.  in  order  to  check  the 
accumulation  of  round-otl  errors  in  the  integration  routine  Since  one  is  free  to  refine  the 
mesh  si/c.  the  trajectory  may  be  made  to  reproduce  itselt  on  each  limp  to  any  accuracy 
(subject  to  machine  constraints)  One  should  use  as  large  a  step  si/c  as  is  compatible  with 
the  accuracy  desired  in  order  to  minimize  computations 

III  INSTRUCTIONS  FOR  USING  RAY  I 

Here  we  describe  the  required  input  data,  output  format,  and  general  interpretation 
of  results  for  RAN  I 

A  INPUT  REQUIREMENTS 

Input  data  is  read  from  cards  under  the  Fortran  free* field  format  This  lonnat  spe¬ 
cification  requires  commas  between  separate  data  items  but  allows  the  items  to  be  positioned 
arbitrarily  on  a  data  card  It  a  single  read  statement  causes  the  input  of  data  continued  on  to 
subsequent  cards,  no  comma  is  required  alter  the  last  item  on  each  card  In  this  case,  the 
card  end  itselt  serves  as  a  legal  delimiter  between  successive  items 

On  the  First  card,  the  integer  variable  ICASI  must  be  punched  This  should  be  tin- 
number  ol  rays  to  be  traced 

On  the  second  card,  and  the  remaining  K'ASI  I  cards,  the  data  items  as  shown  in 
table  2  must  be  punched  Typically  one  card  is  adequate  to  enter  all  items  (or  a  ray. 

B  OUTFIT 

RAY  I  allows  a  quick  check  of  the  input  data,  by  printing  the  variables  for  each  case 
teach  ray)  considered  in  the  form 

(  ASF  I 

OX  *  I  2-02  MFTFRS 

\MA\  *  I  IHOOtXFKH  M 

GRAD  I  -  -  f>J5IO-02  (CHANNFI  WIDTH  GRADIENT) 

(.RAD  2  «  0  0000  H/ (AXIAL  SPFI  D GRADIENT  IN  RANGE) 

CO*  I  0000+00  M/SFC 

AO  *  I  0000+00  M 

ZDOT  *  5  6.^7  -  oi  RADIANS  or  2  0000+01  DEGREES 

ZO  *  0  0000  M 

PROF 'll  I  I  <  l-HIRSCII.  2-PARABOLIC.  3-BIUNEAR) 

It  the  mesh  was  refined  during  any  cycle,  the  next  linets)  will  present  this  information 
in  the  form 

MESH  »  2  AT  STEP  814  . 
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Table  2. 


Card 

Fortran 

Variable 

Name 

Definition  in  RAY  1 

1 

K  ASF 

Nunibei  of  rays  to  be  traced 

‘y 

DX 

Range  increment  for  Kungckutli  raytiaang  routine 
(meters) 

\MAX 

Cut-off  range  tor  riytiacing  (meters! 

(.RADI 

Gtadirnt  in  vet  Ileal  wale  length  4x) 

grad: 

Gradient  in  axial  sound  speed  cglx  1(11/1 

CO 

Initial  axial  sound  speed  (m  sec  I 

AO 

Initial  vertical  scale  factor  (meters! 

/DOT 

Initial  ray  slope  (tanb^l 

ZO 

Initial  ray  depth  (meters) 

IPRINT 

Control  variable  lor  printing  a  detailed  path 
(0  no  .  1  yes) 

ITYPI 

Control  variable  for  SVP  type 

(1  Ihtsch  2  parabolic  .  3  bilinear  1 

3 

DX.  ITYPI 

Data  for  next  tiy 

indicating  one-half  Ihc  input  step  si/e  was  used  past  step  SI4  These  lines  will  not  appear 
if  no  step  si/e  adjustment  was  made 

R  AN  I  employs  a  cut-off  not  only  in  range,  but  also  in  the  number  of  steps  calculated 
(presently  *>999  steps  are  done l  and  in  the  channel  width  parameter  The  user  may  mad- 
vertenlly  set  up  a  constricting  sound  channel  such  that  aixi  becomes  zero  and  then  negative 
in  the  allowed  range  Ray  tracing  is  meaningless  for  such  a  profile  Rather  than  "bombing'' 
under  sush  conditions.  RAN' I  will  test  aixl.  find  it  is  nonpositive,  print  a  message  of  the  form 

CHANNFI  WIDTH  1 1  SS  THAN  0  01  AT  ST  I  I*  704  . 
terminate  ray  (rasing  of  the  present  ra> .  and  go  on  to  the  next  case,  il  there  is  one 

After  RAY  I  has  completed  tracing  a  given  ray  to  YMAX  (or  tor  W)  steps,  or  until 
.Kx>*»0.  whichever  ossurs  first!,  it  prints  the  summary,  beginning  with  the  number  ol  axis 
srossings  (In  contrast  to  RAY2  where  the  ocean  surface  is  at  r  *  0.  RAY  I  always  sets  r  *  0 
at  the  channel  axis  )  Then,  for  each  axis  crossing,  the  following  data  arc  pnntcd 

STI  P  RANOP  TIMI  AVGfDX/DT)  A/AO  MAX  Z  XP'AO  J 

NO  (M>  (SFO  (M/SEC)  <M)  (SEC) 

Here  AO  denotes  the  initial  value  of  the  channel  width  parameter,  and  A/AO  the  dimension¬ 
less  local  value  of  at  x>  This  is  the  same  function  a(  xl  appearing  in  cq  ( 1  2)  above.  For 
constricting  channels,  the  numbers  in  the  A  AO  column  should  be  monotomcallv  decreasing. 


In 


It  is  worthwhile  noting  that  the  AVGtDX  DT I  column  is  derived  from  pairs  of  axis 
crossings  and  the  first  value  corresponds  to  the  average  range  rate  or  horizontal  velocity 
between  the  initial  ray  position  and  the  first  turning  point  For  this  reason,  it  may  not  represent 
an  average  between  turning  points.  For  the  same  reason,  the  last  value  is  arbitrarily  set  to 
zero,  there  being  no  point  beyond  the  last  turning  point  with  which  to  calculate  c.  Similar 
comments  apply  to  the  entries  in  the  J  column,  the  (dimensionless)  cycle  penod  column 
\P  AO.  and  the  maximum  depth  column 

IV  INSTRUITIONS  FOR  USING  RAY2 

In  this  section  we  present  information  required  to  understand  and  run  RAY'2. 
beginning  with  a  brief  description  of  the  necessary  input  parameters  and  the  appropriate 
format  for  these,  the  output,  and  details  of  the  program’s  s«.ope  and  internal  logical 
structure 

A  INPUT  PARAMETER  LIST 

RAY2  consists  of  routines  tor  tracing  rays  in  user-provided  SVPs  A  cubic  spline*' 
is  used  to  smooth  the  input  data  and  to  calculate  interpolated  values  and  derivatives  of  c(z) 
required  by  the  Rungr-kutta  routine  Output  is  provided  in  a  summary  form  similar  to  that 
ot  RAN  I  IMailed  instructions  for  using  RAY'2  are  provided  in  section  IV.  and  section  V 
contains  a  listing  on  the  Fortran  source  code  All  input  data  must  be  read  Iroin  cards  Since 
treetield  format  is  used  tor  input  (with  one  exception,  an  alphanumeric  header  card)  on  may 
simply  place  data  items,  in  sequence,  in  any  fields  on  the  cards  Commas  must  appear 
between  successive  items,  except  when  the  next  card  is  a  continuation  of  the  same  read  list 
items  Ibis  often  occurs  lor  the  SVP  data  In  this  vase,  the  end  of  the  card  itself  avts  as  3 
legal  delimiter  between  the  last  item  on  the  given  card  and  the  first  item  on  the  next  card, 
and  no  comma  is  rrquirrd 

KAY  2  requires  the  input  shown  in  table  3  immediately  alter  the  usual  ,aXQT 
KAY'2.\BS  statement 

Idle  variable  NOSYPS  tells  RAY'2  how  many  SVP  functions  are  to  follow  If  this 
vard  is  incorrectly  punched  tor  misread).  RAY'2  will  either  end  its  computations  before 
rravhing  the  last  data  set  or  attempt  to  read  beyond  the  last  data  set  provided.  If  this  card 
is  missing.  RAY'2  will  recognize  NPTS  as  NOSYPS.  and  then  fail  (probably  with  a  fatal 
execution  erTorl  when  attempting  to  read  vard  .i  under  the  format  2IIO,2Ab 

On  vard  2.  NPT S  is  the  number  of  pairs  of  depth  (ZSVPl  and  sound  speed  (CSVPt 
values  comprising  the  first  SYP  It  NPTS  is  ttiispunvhrd  RAY  2  may  attempt  to  read  past 
the  last  card,  or  may  stop  before  reading  the  last  card,  and  proceed  to  recognize  card  j*l 
as  SYP  data  or  SVP  data  as  card  )♦  I  Results  of  such  input  errors  could  include  many  forms 
of  abnormal  termination 

Immediately  following  NPTS  is  U  NITS,  which  indicates  the  choice  of  units  made 
by  the  user  for  SYP  data  Set  IUNITS  equal  to  1  if  the  units  are  to  be  meters  and  meters 
per  second  and  2  if  feet  and  feet  per  second 

h.  Mathpav.  Innav  Corp  (SPI  SI  and  SP1N2  are  documented  in  a  number  of  manual*  covenng  rendent 
subroutines  I 


I  able  3. 


Card 

Number 

Fortran  Variable  or 
Arias  Name 

Input  Fonnai  Sp<\ 
(Variable  Type) 

1 

NOSVPS 

Free  Field  llntegei) 

T 

N'PTS,  Il'NITS.  MSC.tll 
MSGO! 

:il0.2An 

i 

ZSVPt  1 1.  CSVPt  1 1 
ZSVPtli.  CSVPt  it 

Free  Field  iReal) 

4 

ZSVPt !♦  Ik  CSVPt  1*1!, 

F  rec  Field  IReal) 

J 

ZSVPt  SPIN! 

( SNPlNPTSl 

Free  1  leld  iReal! 

l>\  ZO.  /MAX.  IPLOT 

Free  Field  iReal.  Real. 
Real.  Inlegei  I 

data  for  next  SVP. 
bcpnnini;  with  N'PTS 

(UNITS,  mm.iI). 

MS(,<2(.  ami  ending 

k 

with  l>\.  ZO.  cu 

k+l. 

rial  a  for  Iasi  of 

etc. 

NOSVPS  data  vetv 

MS(>i  1 1  and  MSfi<2i  provide  j  title  lor  the  SVP  to  follow,  consisting  of  12  characters 
selected  from  the  full  range  of  l  MV  A<  I  1 10  I  ortran  and  entered  in  columns  21-32  of  card 

Note  that  card  2  is  read  under  a  tormatted  input  specification. 

Beginning  with  card  3.  and  continuing  on  subsequent  cards  i(  necessary,  is  the  Tirst 
set  of  SVP  data.  As  many  pairs  may  be  entered  on  one  card  as  space  permits  and  pairs  of 
tZSVP.  CSVPi  values  may  be  split,  with  /SVPtii  ending  one  vatd  and  CSVPt  n  beginning 
the  next  card. 

Four  variables  which  control  ravtracing  arc  entered  on  the  card  after  the  last  SVP 
data  card.  l)X.  ZO.  and  ZMAX  arc  the  range  increment  for  the  Rungc-Kutta  algorithm,  the 
initial  ray  depth,  and  the  maximum  depth  allowed  for  a  ray.  respectively  These  variables 
are  interpreted  hv  RA  Y2  as  being  expressed  in  meters  IPLOT  should  be  set  to  0  if  one 
desires  printed  output,  and  I  if  plots  arc  desired  in  addition.  Typical  values  for  L)X  should 
be  approximately  one-sixteenth  of  a  ray  cycle  period,  for  coarse  results,  and  less  for  more 
accurate  raytracing.  If  a  ray  exceeds  depth  ZMAX.  the  remaining  (steeper!  rays  arc  not 
traced  and  RAY2  attempts  to  read  the  next  SVP  data  set 

B  OUTPUT 

RAY2  is  designed  to  allow  the  user  an  immediate  check  of  all  input  The  first  line 
printed  after  the  '••XQT  statement  should  agree  with  the  input  variable  NOSVPS.  and  reads 

IS 


•••  I  HIS  RUN  CONSISTS  OF  tNOSVPS)  SVPtSl  ••• 


V 


ITie  jocund  output  line  consuls  of  the  I  2  character  message  vtnng  MS(.<  1 1.  MSG<2) 
imbedded  as  follows 

•••  SVP  <MS(,il >.  MSG<2)»  ••• 

Dus  simply  identifies  the  particular  SVP  being  computed  It  several  SVPs  are  to  be  run.  one 
such  header  message  appears  at  the  start  of  the  output  tor  ravh  SVP 

The  header  message  is  followed  immediately  h\  a  line  that  vhcvkx  the  vhoice  of  units 
for  the  SVP  data,  and  reads 

USFR  INPIT  UNITS  FOR  THF  SVP  1 1  or  :t  ( I  M  M  SFC  :  FT.  FT  SFO 

UNITS  and  the  units  choice  <1  or  2)  should  agree  Next  the  SVP  is  printed  as  a  pair  of 
columns  labeled  1)1  PTH  and  c  There  should  be  NP1S  pairs  ot  salues 

After  the  SVP  data.  RAY2  prints  the  control  parameters  specifying  how  rays  are  to 
be  traced,  namely  l)X.  /(>.  and  /MAX 

Beginning  with  the  line 

(  ASF  NUMBFR  I  . 

RAY2  prints,  as  presently  structured.  40  rav  patlr  summaries,  for  the  rays  which  are  launched 
Irom  depth  /()  at  the  angles  -0.5°.  -1.0  .  -1.5  .  .  -20  .  with  respect  to  horizontal  F.ach 

vase  is  summarized  by  a  set  of  initial  conditions  printed  in  the  form 

(  ASF  NIMBI  R  1 

INITlAl  ANGI.F  *  0  5  l)F(.RI  I  S 

I N  IT  I A I  SOUND  SPHFI)  *  1495  d  M  SFC 

PRI  OR  FI  I)  VI- RTF  XING  SPI  10=  1  <00  2  M  SI  C 

The  first  and  second  lines  are  redundant  since  the  initial  angle  associated  with  a  given  case  is 
alwass  one-halt  the  vase  number  Nevertheless,  it  is  useful  to  sec  at  a  glance  the  angle  of  the 
ray  being  considered,  and  if  internal  vhanges  are  made  to  produce  a  different  angular  interval 
for  rays,  these  two  lines  would  sene  as  checks  on  the  code  changes  and  would  simplify  some 
data  analy  sis 

The  initial  sound  speed  is  computed  by  the  spline  curse  fitting  routine  SPLNI.h  and 
used  to  compute  a  predicted  vertex  speed  according  to  Snell  s  law. 

‘vertex  *  rt*0,/corf  • 

F.ach  case  is  also  summarized  by  a  table  showing  the  history  of  ray  turning  points. 

As  presently  written  RAY2  automatically  traces  rays  through  four  turning  points.  The 
variables  printed  at  each  turning  point  (with  the  exception  of  J  for  the  first  turning  point 
and  the  average  speed  for  the  last  turning  point)  are 

K.  XTP.  TTP.  ZTP.  (TP.  J(K).  AVGSP  . 
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Diese  stand  for  the  turning  point  number,  k.  the  range  to  the  kth  turning  point,  the  time 
elapsed  at  the  Kth  turning  point,  the  depth  of  the  Kth  turning  point,  the  ray  speed  at  the 
Kth  turning  point,  the  action  integral  /  sind  dt/c  computed  over  the  last  one-half  cycle,  and 
the  average  horizontal  ray  speed  from  the  Kth  to  the  t  K ♦  I  Kt  turning  point 

f  igure  5  illustrates  the  physical  significance  of  these  summary  data  items  Note 
that  the  first  value  printed  in  the  J  column  is  really  an  incomplete  integral,  evaluated  from 
the  starting  depth  zq  to  the  first  turning  point,  and  should  he  ignored  KAY2  computes  the 
average  horizontal  ray  velocity  between  pairs  ol  turning  points  The  first  AVGSPvaluc 
pnnted  out  is  assixiated  with  the  first  and  second  turning  points,  and  the  fourth  value  is 
automatically  set  to  zero 


figure  S  A  sketch  of  lire  phcucal  meaning  of  the  summary  data  items  in  KAY2 

V  INTERNAL  LOGICAL  STRUCTURE  OF  RAY  I  AND  RAY2 

In  this  section  we  shall  describe  the  logical  structure  of  the  Fortran  routines 
comprising  RAYI  and  RAY2  Special  emphasis  is  placed  on  the  areas  in  which  slight  modi¬ 
fications  m  these  codes  would  enhance  their  usefulness,  particularly  in  the  direction  of  per¬ 
mitting  better  output  (eg  .  plots)  and  more  flexibility  through  user -selected  input  variables 
bey  ond  those  now  available 

A.  INTERNAL  STRUCTURE  OF  RAYI 

As  it  is  presently  written.  RAYI  contains  comment  cards  defining  most  of  the 
variables  of  interest  (see  section  VI).  A  flowchart  of  RAYI  appears  in  figure  6  1C ASF  is 


HEAD  ICASI 


f  igurc  t>  Howchut  of  RAVI 


the  number  of  eases  to  be  considered  lach  case  consists  of  a  particular  choice  for  an  SVP. 
initial  ra>  depth  and  slope,  and  parameters  for  the  range  step,  maximum  range,  and  gradients 
of  the  range-dependent  parameters  (axial  sound  speed  and  channel  width  parameter  in  cq 
it'll  In  addition.  RAVI  reads  IPRIVT  which  allows  or  suppresses  a  longer  printout  of 
the  trajectory 

An  extremely  valuable  addition  to  RAVI  would  be  a  plotting  routine  which  simply 
plots  the  ray  paths  One  could  then  observe  immediately  (without  calling  the  subroutine 
PA  I II  >  the  effect  of  squeezing  or  expanding  SYPs  on  ray  paths  A  second  much  needed 
improvement  consists  of  permitting  several  rays  to  be  traced  in  a  given  range-dependent  SVP. 
without  the  requirement  of  punching  all  SVP  data  and  other  parameters  anew  for  each  ray. 

The  subroutine  TRAJPC  computes  a  ray  path,  including  the  time  and  average  hori¬ 
zontal  velocity  at  each  turning  point,  and  the  action  integral  over  each  half-cycle.  Beginning 
at  line  'b  (see  section  V|>  TRAJtC  computes  the  new  ray  slope  (DPRIZi.  These  instructions 
end  at  line  40  of  TRAJKC.  and  consist  of  the  Rungc-Kutta  algorithm  mentioned  earlier 


Oner  an  advanced  v aluc  of  /  is  found.  I  RAJfcC  finds  the  new  sound  speed  at  this 
depth  and  range,  which  in  turn  is  used  to  compute  a  time  increment  given  by  At  =  As/c,  and 
a  new  increment  of  the  action  integral  given  by  c“*  suit)  A  /  These  steps  correspond  to 
lines  42  through  55.  I  wo  tests  are  made  here  to  insure  that  the  user  has  not  allowed  the 
channel  width  parameter  to  vanish,  and  to  avoid  any  illegal  tor  rather  invalid)  reference  to 
array  elements  b>  looking  backwards  from  the  very  lust  step  of  the  ray  tracing 

Hie  remaining  portions  of  TR.AJFC  handle  axis  crossings  and  adjust  the  step  si/e 
should  there  be  fewer  than  I  50  steps  per  half-cycle  \5e  recommend  that  an  input  variable 
be  added  to  permit  varying  this  particular  dunce  at  run  time 

In  section  II  we  derived  the  filler- l.agrange  equation,  cq  t1).  for  Fermat's  variational 
principle  RAN  I  calculates  the  right-hand  side  ol  this  equation  ol  motion  m  the  function 
subroutine  F.  which  is  called  by  I  R  \ J I  (  Ihe  three  SVI\  m  equations  (12)  lead,  respec¬ 
tively  .  to  the  results 


.  s  x  s  „  .X 

F  ■  -2(1  ♦  /  *)  /  (a*  ♦  /*  I  ♦  /(I  ♦/‘igx  cq 

•  ,  N  >  1 

-2  /t  I  ♦  r~  a)  gj  (a“  ♦  /- )  , 

I  -  1/  al"  I  I  ♦  /*  )  ( I  -  1/  at-) 

•  ,  1 

♦  /(  1  ♦  /■  )  gs  Cq 

•  ,  N  >  J  ">  ■> 

-  /<  I  ♦  /*  </*  s  )  g|  1  1  -  /*  a"  >  . 

.  y  , 

l  Ml  ♦  /  '  )  I  .1  <  /!)♦/(!  ♦  /  ‘  )  g  s  Cq 

-  /(  I  ♦  /*  i  /  gj  ala  ♦  7ll  .  /  v  0  (  >  0 1  . 


(  I.Ti l 


(13b) 


( 13c ) 


as  one  can  verify  easilv  bv  direc  t  differentiation  Hcrvgj  and  gs  are  respectively  da  dx  and 
dci/=0)  dx 

MAX/  and  PA  I  Mare  short  routines  Ihe  former  searches  each  half  cycle  of  the 
trajectory  and  returns  a  list  ol  the  maximum  off-axis  distance  reached  by  the  rav  Ihe  latter 
pnnts  250  points  equally  spaced  along  the  trajectory,  or  as  many  points  as  arc  found  before 
the  maximum  range  is  reached 

RAVI  should  be  expanded  to  handle  a  wider  class  of  range-dependent  analytic  S\  IV 
and  to  permit  connecting  range-independent  SVPs  smoothly,  as  discussed  in  section  II 


B  INTERNAL  STRIKTIJRI  OF  RAY 2 

RAY 2  consists  of  only  two  routines  and  two  cubic  spline  routines  resident  on  the 
NOSC  UNIVAC  I  1 10. h  Unlike  RAVI.  RAN  2  does  provide  a  plot  subroutine  and  uses  any 
user-provided  SVP  We  recommend  that  the  plot  routine  calls  be  modified  to  allow  several 
rays  to  be  plotted  on  a  single  graph  Also  desirable  would  be  a  plot  of  the  user-input  SVP 
as  smoothed  by  the  spline  routines 

A  flowchart  of  R  AN' 2  is  given  in  figure  7  Although  somewhat  complicated  in 
appearance,  the  flowchart  involves  only  three  basic  loops  The  outermost,  which  corresponds 


1  7 


to  the  DO  loop  with  object  l>0  in  R A V 2  (line  lo.l  i,  steps  through  the  NOSVPS  sound  speed 
profiles  to  be  considered  I  he  next,  with  object  bO  (line  lo2>  controls  the  40  ravs  to  be 
traced  in  each  SVP  Ihe  innermost  loop  checks  the  number  of  lull-cycles  traced  and  jumps 
to  a  new  initial  ra>  angle  when  this  exceeds  or  equals  live.  Hence  each  ray  is  traced  over 
fixe  turning  |>omls  \s  m  RAN  I,  a  number  of  tests  are  made  at  each  step  to  pick  up  the 
turning  points,  except,  ot  course,  that  surlace  relies  lions  as  well  as  lotal  internal  reflections, 
may  occur  here 

One  should  note  that  it  the  depth  limit  is  reached  In  an>  ray.  RAN  2  will  automatically 
print  a  message  to  this  effect  and  look  lor  a  new  S\  I*  It  will  not  trasc  the  next  ra\  (0  5 
degrees  steeper  l  in  the  present  SV|* 

RAN  2  would  benelil  (roni  a  revision  allowing  a  variable  spacing  between  rays  and 
a  variable  number  ol  ravs,  instead  ol  0  ^  and  40.  respectively 


VI  I  1ST IM. SOI  RAVI  ANDKAV2 

In  this  section  we  present  listings  (computer  printout  I  ol  the  I  ortran  source  code 
comprising  RAN  I  and  RAN  2  Ihe  reader  vhoulsl  note  that  in  some  vases  pagination  becomes 
awkward  when  converting  printout  to  report  format  In  particular,  separate  pages  vannot 
always  reler  to  separate  routines  Ihe  I  ortran  line  sequence  numbers  which  appear  in  the 
lettdiand  column  should  therefore  be  used  to  follow  a  routine  acri*ss  several  pages  or  to 
distinguish  multiple  routines  from  eaJi  other  on  one  page  In  addition,  the  line  sequence 
numbers  are  used  elsewhere  in  this  report  to  refer  to  particular  sections  ol  the  codes 
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